Identification of drought-tolerant hub genes in Iranian KC-2226 genotype of Aegilops tauschii using transcriptomic analysis

Aegilops tauschii, as a donor of D genome to the bread wheat with a valuable source of resistance to different biotic and abiotic stresses, is used to improve the quality of wheat cultivars. Every genotype has a specific genetic content, the investigation of which can lead to the identification of useful genes such as stress tolerance genes, including drought. Therefore, 23 genotypes of Ae. tauschii were selected to evaluate their morphological and physiological traits under greenhouse conditions. Among them, a superior tolerant genotype (KC-2226) was chosen for transcriptomic analysis. Our result showed that 5007 and 3489 genes were deferentially up- and downregulated, respectively. Upregulated genes were involved in photosynthesis, glycolysis/gluconeogenesis, and amino acid biosynthesis whereas downregulated genes were often engaged in DNA synthesis, replication, repair and topological changes. The result of protein–protein interaction network analysis showed that AT1G76550 (1.46), AT1G20950 (1.42), IAR4 (1.19), and PYD2 (1.16) among upregulated genes and THY-1 (44), PCNA1 (41) and TOPII (22) among down-regulated genes had the highest interactions with other genes. In conclusion, Ae. tauschii employs elevated transcription of specific genes involved in photosynthesis, glycolysis and gluconeogenesis and amino acid biosynthesis pathways rather than genes active in DNA synthesis and repair to provide the energy needed for the plant to survive under stress conditions.

RNA extraction and sequencing. Total RNA was extracted from 200 mg ground leaves samples (4 leaf stage) using RNeasy Plant Mini Kit (Qiagen, Germany) according to the manufacturer's protocol. The quantity and quality of the extracted RNA were analyzed using NanoDrop 1000 spectrophotometer (Thermo Scientific, USA) and 1% agarose gel electrophoresis, respectively. The extracted RNA samples were treated with DNase I (Thermo Scientific, USA) to eliminate probable genomic contamination. Each of the control and drought stress treatments had three replications. Subsequent quality control of the extracted RNA was done using a QC Bioanalyzer (Agilent Technologies, CA, USA). The Poly-A selection, preparation of cDNA, adapter ligation, clusters formation, and sequencing was carried out at Beijing Genomes Institute following the manufacturer's instruction, using TruSeq Stranded total RNA with Ribo-Zero Plant kit (Illumina, USA). The sequencing was performed by Poly-A mRNA Capture method using a Nova-seq 6000 platform to produce 6 GB raw data of 100 bp paired-end reads.
Morphological and physiological analysis of genotype KC-2226. According to Table 3, The RWC, seedling root length, and percentage of yellow leaves were decreased whereas seedling shoot length, the number of leaves and tillers, fresh and dry weight of roots and shoots as well as the percentage of tubular leaves were increased in the seedlings of genotype KC-2226 under drought stress. In mature plants, shoot dry weight and spikes number increased under drought stress, but root dry weight, and the length of flag leaf, peduncle, and spikes decreased (Table 3). In addition, under drought stress, proline content and activities of SOD, CAT, APX, and POX enzymes increased (Table 3).
Data quality and mapping. A total of 111,181,350 raw reads were constructed from two cDNA libraries.
After trimming and removing the low-quality reads, a total of 95,042,354 clean reads were obtained (Table 4). A total of 89,089,860 reads were mapped to the genome with both forward and reverse primers, which included 93.75% of the total clean reads. A total of 3,792,499 reads (3.96% of clean reads) were mapped to genome with one either forward or reverse primer, and finally, a total of 2,159,995 reads (2.28% of clean reads) were not mapped to the genome (Table 4).   2). According to the volcano plot, the distribution of upregulated and downregulated genes is almost similar. In this figure, the genes that were significant at the 0.05 probability level are shown. In total, 5007 genes were upregulated whereas 3489 genes showed a downregulation under drought stress (Fig. 2).    (Figs. 3 and 4). Upregulated genes were often involved in biological activities of cellular component disassembly, carboxylic acid and organic acid catabolic process, cellular amino acid biosynthetic and metabolic process with enrichment values from 2.91 to 10.01, which were active in chloroplast, cytosol, and membrane. The molecular functions of these genes were mainly catalytic, lyase, and ion binding with enrichment values from 1.4 (cellular anatomical entity) to 18.22 (for oxidoreductase activity) (Fig. 3). Downregulated genes were involved in biological processes of deoxyribonucleotide biosynthetic process, DNA unwinding cell cycle checkpoint signaling, and DNA conformation change with the enrichment values from 11.79 to 25.6. The deoxyribonucleotide biosynthetic process and DNA unwinding involved in DNA replication had the highest enrichment values (25.6 and 23.72, respectively).
These genes were active in nuclear replisomes, nuclear replication fork, DNA polymerase complex and chromosomes, which were involved in molecular functions including DNA-methyltransferase activity, transporter activity, DNA helicase activity and ATP-dependent activity with the enrichment values from 5.24 to 31.29. The DNA-methyltransferase activity had the highest enrichment value (31.29) in all gene ontology components (Fig. 4).
Our result showed that MAPK signaling pathway is one of important active pathways involved in drought stress conditions (Fig. 6). Once stress signals were perceived by cell membranes, this pathway is regulated by PRY/PRL as soluble abscisic acid (ABA) receptors in an ABA-dependent way. The 2C-type protein phosphatases (PP2C), which inactivates SNF1-related protein kinases2 (SnRK2) through dephosphorylation, is inhibited in the presence of ABA.    Table 5. The number of nodes and edges were approximately two and eight times, respectively, more in PPI network of downregulated genes compared with those in upregulated genes. The average local clustering coefficient of PPI network of up-and downregulated genes were 0.332 and 0.41, respectively, with the same P-value. Gene network was constructed through the data including nodes and edges acquired from PPI analysis (Fig. 7). Genes pyrophosphate-fructose 6-phosphate 1-phosphotransferase subunit alpha: PFP-ALPHA (1.46),    (18), THY-1 (16) and RNR1 (15) among downregulated genes had the highest scores and recognized as hub genes (Fig. 7).
Expression analysis and validation of RNA-seq data. The expression level of two upregulated genes IAR4 and PYD2, and two downregulated genes TOPII and THY-1 under control and drought stress treatments were analyzed and then compared with RNA-seq results. The expression of IAR4 and PYD2 genes increased under the control and stress treatments with relatively higher expression in the drought stress treatment than those in the control. The expression of IAR4 (6.28) and PYD2 (3.16) under stress treatment were almost three In gene networks, the color intensity determines the rank of the genes, the higher the color intensity and closer to red, the higher the score, and the closer to yellow, the lower the score. Genes with higher scores are more important to us and are our hub genes. www.nature.com/scientificreports/ times more than their control (2.2 and 1.3, respectively). The expression of TOPII and THY-1 genes decreased under control conditions and under drought stress. The severity of the decrease in the expression of these two genes in drought stress conditions was significantly higher than in the control. The expression of TOPII gene under drought stress decreased by 4.5 times compared to the control condition, while the expression of THY-1 gene decreased by 5.6 times under drought stress compared to the control condition (Fig. 8). The RT-qPCR result for selected genes was almost consistent with the RNA-seq data. LogFC of IAR4, PYD2, TOPII and THY-1 genes were 13.94, 4.54, −4.2 and −3.58, respectively.

Discussion
Since the DD genome of Ae. tauschii Coss is a good source of resistance to different biotic and abiotic stresses, resistance genes from this species can be transferred to bread wheat 16 or its close relatives through classical approaches for breeding purposes and increasing the quality of wheat cultivars 5 . To find the most important genes involved in drought stress, at first, genotype KC-2226 was selected among 23 other Iranian Ae. tauschii genotypes as the most tolerant one according to morphological and physiological properties. Transcriptomic analysis of genotype KC-2226 was carried out and DEGs were identified under drought stress. PPI analysis was done and the gene network was constructed using CytoHubba by employing MCC method, which had a better performance on the prediction accuracy of important proteins from the PPI network. The most important upregulated genes were AT1G76550, AT1G20950, IAR4, and PYD2 whereas the most important downregulated genes included THY-2, PCNA1, and TOPII (Fig. 7).
Genes AT1G76550 and AT1G20950 are involved in the process of photosynthesis and response to glucose. These genes are located in cytoplasm and mitochondria and are a part of pyrophosphate-dependent phosphofructokinase complex and alpha-subunit complex. These genes are active in the glycolysis/gluconeogenesis, fructose and mannose metabolism and pentose phosphate pathways 33 , having molecular functions including ATP binding, phosphofructokinase, and diphosphate-fructose-6-phosphate 1-phosphotransferase activity 34 .
The PFPs are involved in response to different stresses including dehydration, high salt, phosphate starvation anoxia and wounding 35,36 . They play roles in gluconeogenesis, glycolysis, stabilization of triose-phosphate and hexose-phosphate pools, and regulation of inorganic pyrophosphate (PPi) concentration during synthesis and degradation of sucrose as well as adaptability to stresses 36-38 . Lim et al. (2014) investigated the response of PFP double and quadruple knockout mutants to osmotic and salt stresses to clarify the role of PFP in the stress tolerance of Arabidopsis seedlings. The expression of PFP subunit genes increased in response to salt and osmotic stresses. These findings suggest that PFP plays a role in adapting to salt and osmotic stresses 39 .
The IAR4 is involved in the biosynthesis of acetyl-CoA from pyruvate, glycolytic processes as well as auxin conjugate sensitivity and homeostasis in root development 40 . This gene is located in the cytosol, mitochondria, and mitochondrial matrix with pyruvate dehydrogenase (acetyl-transferring), and cobalt and zinc ion binding activities 41 . It is engaged in glycolysis/gluconeogenesis, pyruvate and carbon metabolism and citrate cycle pathways 33 . Fu et al. (2019) reported that IAA-CONJUGATE-RESISTANT 4 (IAR4) plays a key role in primary root growth under salt stress conditions. Mutation of IAR4 led to increased sensitivity to salt stress conditions, with strongly inhibited primary root growth and reduced survival rate in two iar4 mutant alleles 42 .
The PYD2 is involved in the beta-alanine biosynthetic process, cellular response to nitrogen levels, pyrimidine nucleobase catabolic process, and uracil catabolic process. The PYD2 is located in endomembrane system, endoplasmic reticulum, Golgi apparatus, and plastid with dihydropyrimidinase activity and metal ion binding activities. It is engaged in amino acid and beta-alanine biosynthesis pathways.
The THY-1 is involved in the biosynthesis of 10-formyltetrahydrofolate and dTMP, methylation and onecarbon metabolic process. This gene is in cytosol and mitochondria with dihydrofolate reductase and thymidylate synthase activities, which are crucial for DNA synthesis. These two enzymatic activities in plants are expressed as one bifunctional enzyme 43 . In addition, it is active in pathways of cofactor biosynthesis and tetrahydrofolate biosynthesis. Gorelova et al. (2017) showed that one of the DHFR-TS (dihydrofolate reductase-thymidylate synthase) isoforms (DHFR-TS3) operates as an inhibitor of its two homologs, thus regulating DHFR and TS activities and, as a consequence, folate abundance. In addition, a novel function of folate metabolism in plants  www.nature.com/scientificreports/ is proposed, i.e., maintenance of the redox balance by contributing to NADPH production through the reaction catalyzed by methylenetetrahydrofolate dehydrogenase, thus allowing plants to cope with oxidative stress 44 . The PCNA1 is involved in the biological elongation of the leading strand, repair of mismatch, regulation of DNA replication, regulation of cell cycle and translation. It plays different functions including DNA binding, DNA polymerase processivity factor activity in the cytoplasm, cytosol, nucleolus, nucleus and PCNA complex. The PCNA1 is active in DNA replication, base and nucleotide excision and mismatch repair, cell cycle 45 . Ghabooli et al. (2013) conducted a proteomics study to understand the molecular mechanisms underlying water stress tolerance induced by Piriformospora indica in barley. They reported that the abundance of PCNA decreased in response to drought stress. However, P. indica colonization resulted in an increase in the abundance of this protein under drought conditions 46 .
The TOPII is engaged in DNA topological change and as an active cellular component of intracellular membrane-bounded organelle has different roles including ATP binding and hydrolyzing, metal ion binding, DNA binding, and double-strand break. Topoisomerases mitigate topological stress by untangling and relaxing the supercoiled DNA in both eukaryotes and prokaryotes 47,48 . John et al. (2016) over-expressed topoisomerase II (TopoII) in tobacco (Nicotiana tabaccum) and examined its role in growth and development as well as salt (NaCl) stress tolerance. They revealed that NtTopoII1-α over-expression in tobacco confers salt stress tolerance to the transformed lines as compared to wild-type plants. TopoII over-expression changed the morphology of the transgenic plants and improved the seed germination on a salt-supplemented medium 49 .
Altogether, upregulated genes were often involved in photosynthesis, glycolysis and gluconeogenesis, amino acid and beta-alanine biosynthesis whereas downregulated genes were mainly involved in DNA synthesis, replication, repair and topological changes. These results indicated that photosynthesis, glycolysis, and gluconeogenesis are prioritized in Ae. tauschii grown under drought stress to provide more energy for the plant than DNA-related activities. Still, understanding the mechanism of the most important pathways dealing with drought stress needs further studies.

Conclusion
Ae. tauschii is known as a plant tolerant to all kinds of biotic and abiotic stresses. Therefore, identifying the genes with the highest interaction, both up-and downregulated genes, under drought stress is of prime importance. According to our results, the plant exploits transcription of specific genes (photosynthesis, glycolysis and gluconeogenesis, amino acid, and beta-alanine biosynthesis pathways) instead of DNA synthesis and repair under stress conditions to provide the energy needed for the plant to survive under stressful conditions.

Data availability
The datasets generated during the current study are available in the NCBI database, BioProject ID of PRJNA868361.